
inde_dilution  = 1;
l_reallocation = 1;
demand_ext     = 1;
shiftdist      = 1;

k0_vec   = [0,sub_time];
incid_vec= [0.1:0.1:1];
k_vec   = (1-sub_size)*ones(1,length(incid_vec));
ip2_sub=100;ip1_sub=100;
W_mat=zeros(length(k0_vec),length(k_vec),1);

str_region="South_";
str1="params_";
str_3="calibration_Table2";
str_save=append(path_calibration,str1,str_region,str_3,".mat");
filename=cell2mat(str_save);
load(filename)
read_pars;
vars=zeros(3,1);
vars(1)=pars(30);
vars(2)=A_s_ben;
vars(3)=pars(30);
n_p=length(pars);
pars(n_p+1)=incid; pars(n_p+2)=k;
fun_FEFC;

varzeta_s_ben=(r+delta)*Ev0;
k_s_ben0=varzeta_s_ben/(r+delta)+ B0*Ex0;
pars(8)=varzeta_s_ben;
pars(n_p+2)=k_s_ben0;k_s_ben=k_s_ben0;
A=A_s_ben;read_pars;read_pars_s;
Z_s=z_0_s*reshape(exp(z_vec)'*ones(1,N_b),Ntot,1);
options = optimoptions('fsolve','OptimalityTolerance',1e-10);

%% compute nwe steady state and dynamics
welf_disc_all=[];
ip2=0;inde_solved = ones(length(k0_vec),length(k_vec));
while ip2<length(k_vec)
    ip2
    ip2=ip2+1;
    size_subsidy=1-k_vec(ip2);
    incid = incid_vec(ip2);
    ip=0;
    while ip<length(k0_vec)
        ip=ip+1;
        ip
        global cub_spl_v cub_spl_x pp_star Amin

        %find the equilibrium profitability of firms that is consistent withthe entry cost
        k_s = (1-size_subsidy*(1-k0_vec(ip)))*k_s_ben0;
        k0_s = k0_vec(ip)*size_subsidy*k_s_ben0;%part ex-post
        pars(7)=k_s;pars(10)=k0_s;pars(n_p+1)=incid;pars_s=pars;

        %% find new value added and debt at entry
        n_try=0;
        fv=ones(3,1);A_0=A_s_ben*0.99; if k0_vec(ip)>0.4,B_0=1.2*B0_s;else,B_0=0.25*B0_s;end
        if ip==1,vars_0=[B_0,A_0,B_0];end
        while max(abs(fv))>1/10000 && n_try<20
            n_try=n_try+1;
            [vars,fv]=fsolve('fun_EntryEq',vars_0,options,pars_s);
            if incid==1
                A_0=A_0*0.95;B_0=B0_s;
                vars_0=[B_0,A_0,B_0];
            else
                vars_0=vars;
            end
        end

        B0_s     = max(0,vars(1));A_s_new=vars(2);B00_s=max(0,vars(3));
        if max(abs(fv))>1/1000 ||  B00_s/(z_0*A_s_new)/zlow_s>ppval(pp_bar,A_s_new)
            warning('no solution'),
            inde_solved(ip,ip2)=0;
        end

        b_common=B0_s/(z_0*A_s_new);
        b0_low_s=b_common/zhig_s;
        b0_mid_s=b_common/zmid_s;
        b0_hig_s=b_common/zlow_s;
        %compute debt buyback when rebate is expost
        if k0_s>0 && b0_hig_s>ppval(pp_star,A_s_new)
            b0_hig_s=fun_buyback(b0_hig_s,A_s_new,k0_s/(z_0*A_s_new),cub_spl_v,cub_spl_x);
        end
        if k0_s==0 && size_subsidy==0,b_up_ss=b_up; b_md_ss=b_md;b_dn_ss=b_dn;end
        b_up=b0_hig_s;b_md=b0_mid_s;b_dn=b0_low_s;
        %untreated firms
        b_common0 = B00_s/(z_0*A_s_new);
        b_dn2=b_common0/zhig_s;
        b_md2=b_common0/zmid_s;
        b_up2=b_common0/zlow_s;
        
        
        % find invariant distribution (b,z)
        b_up1=b0_hig_s;b_md1=b0_mid_s;b_dn1=b0_low_s;
        q_hig=q_hig_s;q_low=q_low_s;q_mid=q_mid_s;

        dt = 0.01;
        T  = 5/dt;
        n_sim=25000;

        b_up_fund=b_up_ss;  b_md_fund=b_md_ss;  b_dn_fund=b_dn_ss;
        b_up_unfund=b_up_ss;b_md_unfund=b_md_ss;b_dn_unfund=b_dn_ss;
        b_bar=ppval(pp_bar,A_s_ben);
        get_sim_exit;
        exit_0_t(:,ip,ip2)=exit_prob_yearly;

        b_up_fund=  b_up1;b_md_fund=b_md1;  b_dn_fund=b_dn1;
        b_up_unfund=b_up1;b_md_unfund=b_md1;b_dn_unfund=b_dn1;
        b_bar=ppval(pp_bar,A_s_new);
        get_sim_exit;
        exit_1_t(:,ip,ip2)=exit_prob_yearly;

        b_up_fund=b_up2;  b_md_fund=b_md2;  b_dn_fund=b_dn2;
        b_up_unfund=b_up2;b_md_unfund=b_md2;b_dn_unfund=b_dn2;
        b_bar=ppval(pp_bar,A_s_new);
        get_sim_exit;
        exit_2_t(:,ip,ip2)=exit_prob_yearly;

        b_up_fund=b_up1;  b_md_fund=b_md1;  b_dn_fund=b_dn1;
        b_up_unfund=b_up2;b_md_unfund=b_md2;b_dn_unfund=b_dn2;
        b_bar=ppval(pp_bar,A_s_new);
        get_sim_exit;
        exit_3_t(:,ip,ip2)=exit_prob_yearly;
    end

end

% exit_ss     = mean(mean(squeeze(exit_0_t(1:4,1,:)),1));
% exit_fund   = mean(squeeze(exit_1_t(1:4,1,:)),1);
% exit_unfund = mean(squeeze(exit_2_t(1:4,1,:)),1);
% exit_elig   = mean(squeeze(exit_3_t(1:4,1,:)),1);
% 
% 
% figure(1)
% plot(100*(incid_vec),100*exit_unfund-100*exit_ss,'-b','LineWidth',6),hold on
% plot(100*(incid_vec),100*exit_fund-100*exit_ss,'--r','LineWidth',6),hold on
% plot(100*(incid_vec),100*exit_elig-100*exit_ss,'-.g','LineWidth',6),hold on
% %yline(100*exit_ss,':k','Status quo','LineWidth',3),hold off
% grid on
% yline(0,'--k','Status quo','Interpreter','Latex','Fontsize',18,'LabelHorizontalAlignment','left')
% xline(20,':k','ISS','LineWidth',2),hold off
% %plot(100*(incid_vec),exit_ss*ones(length(incid_vec),1),':k','LineWidth',2),hold off
% xlabel('$\iota$, \%','Fontsize',18,'Interpreter','Latex')
% ylabel('Avg. exit rate years 1-4 in dev. from status quo, \%','Fontsize',18,'Interpreter','Latex')
% le=legend('Unsubsidized','Subsidized','Eligible');
% set(le,'Interpreter','Latex','Fontsize',16,'Location','SouthEast');
% str1="Text/figures/figure_exit_south_tau0.pdf";
% str_save1=append(root,str1);
% str2="Presentation/figures/figure_exit_south_tau0.pdf";
% str_save2=append(root,str2);
% settt
% print(gcf,'-dpdf',str_save1)
% print(gcf,'-dpdf',str_save2)
% 
% 

exit_ss     = mean(mean(squeeze(exit_0_t(1:4,2,:,:)),1));
exit_fund   = mean(squeeze(exit_1_t(1:4,2,:,:)),1);
exit_unfund = mean(squeeze(exit_2_t(1:4,2,:,:)),1);
exit_elig   = mean(squeeze(exit_3_t(1:4,2,:,:)),1);

figure(2)
plot(100*(incid_vec),100*exit_unfund-100*exit_ss,'-b','LineWidth',6),hold on
plot(100*(incid_vec),100*exit_fund-100*exit_ss,'--r','LineWidth',6),hold on
plot(100*(incid_vec),100*exit_elig-100*exit_ss,'-.g','LineWidth',6),hold on
%yline(100*exit_ss,':k','Status quo','LineWidth',3),hold off
grid on
axis([10 90 -4 4])
yline(0,'--k','Status quo','LineWidth',2,'Interpreter','Latex','Fontsize',18,'LabelHorizontalAlignment','right')
xline(20,':k','ISS','LineWidth',2,'FontSize',18,'Interpreter','Latex','LabelVerticalAlignment','bottom')
%plot(100*(incid_vec),exit_ss*ones(length(incid_vec),1),':k','LineWidth',2),hold off
xlabel('$\mathbf{i}$, \%','Fontsize',18,'Interpreter','Latex')
ylabel('Avg. exit rate years 1-4 in dev. from status quo, \%','Fontsize',18,'Interpreter','Latex')
le=legend('Unsubsidized','Subsidized','Eligible');
set(le,'Interpreter','Latex','Fontsize',16,'Location','SouthEast');
str1="/figures/figure_exit_south_tau50.pdf";
str_save1=append(root,str1);
settt
print(gcf,'-dpdf',str_save1)